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ABSTRACT 

We carry out hydrodynamical simulations of galaxy formation that simultaneously 
follow radiative transfer of hydrogen- ionising photons, based on the optically-thin 
variable Eddington tensor approximation as implemented in the GADGET code. We 
consider only star-forming galaxies as sources and examine to what extent they can 
yield a reasonable reionisation history and thermal state of the intergalactic medium at 
redshifts around z ~ 3. This serves as an important benchmark for our self-consistent 
methodology to simulate galaxy formation and reionisation, and for future improve- 
ments through accounting of other sources and other wavelength ranges. We find that 
star formation alone is sufficient for reionising the Universe by redshift z ~ 6. For a 
suitable choice of the escape fraction and the heating efficiency, our models are ap- 
proximately able to account at the same time for the one-point function and the power 
spectrum of the Lyman-a forest. The radiation field has an important impact on the 
star formation rate density in our simulations and significantly lowers the gaseous and 
stellar fractions in low-mass dark matter halos. Our results thus directly demonstrate 
the importance of radiative feedback for galaxy formation. The spatial and tempo- 
ral importance of this effect can be studied accurately with the modelling technique 
explored here, allowing more faithful simulations of galaxy formation. 
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1 INTRODUCTION 



!• ■ In the standard cosmological model, the primordial gas re- 
. combines and becomes neutral around redshift z c± 1000. 
n However, the absence of Gunn- Peterson t roughs in the spec- 

^ ■ tra of high red shift quasars up to z < 6 ([White et al ] l2003l : 
iFan et alJlioO^ ) suggests that hydrogen is highly ionised at 
low redshift. Thus, there must be a period in the history 
of the Universe when hydrogen became ionised again, but 
it is still an open question when the process of this cosmic 
reionisation started, how it proceeded in detail, and which 
sources of radiation were primarily responsible for it. 

An important observational clue about reionisation is 
provided by the total electron-scattering optical depth to the 
last scattering surface of the cosmic microwave background 
(CMB), found to be r es 0.08785 ± 0.0 0072 by the latest 
WMAP7 data release ([Larson et alj |2010). This points to an 
early start and possibly extended period of reionisation. Fur- 
ther observational information about the history of hydrogen 
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reionisation can be inferred through different astrophysical 
phenomena. The Gunn-Peterson troughs in quasar spectra 
are sensitive to small trace amounts of neutral hydrogen dur- 
ing the late stages of reionisation (z ~ 6); Lyman-a emit- 
ting galaxies probe the intermediate stages of reionisation 
(7 < z < 15); the 21-cm line background and gamma ray 
bursts (GRBs) probe the early stages of reionisation, when 
the Universe was mostly neutral (10 < z < 30); and finally, 
the cosmic microwave background (CMB) polarisation pro- 
vides important data on the free electr on column density in - 
tegrated over a large range of redshifts ([Alvarez et alj|2006l ) . 
In the future, upcoming observations from new radio tele- 
scopes such as LOFAR promise to be able to map out the 
epoch of cosmic reionisation in unprecedented detail, mak- 
ing theoretical studies of this process especially timely. 

Population III stars are commonly believed to be the 
first significant sources of photons for reionisation. They 
are very massive, short -lived stars and form around redshift 
z ~ 30, or even earlier (|Gao et 111 120071 ). Another, in princi- 
ple plausible c andi date for causing reionisation, are quasars. 
iMadau et all (| 19991 ) compare luminosity and spacial density 
functions of quasars and star-forming galaxies and show that 
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quasars alone are unlikely to be the dominant sources of pho- 
tons for reionising the Unive r se. On the other ha nd other 
groups (jTrac fe Gnedinl 120091 : iThomas et all l2Q09h predict 
in recent work that quasars can produce enough ionising 
photons to be able to reionise the Universe. 

In many scenarios, reionisation is assumed to be driven 
mainly by high- m ass objects with mass M > 1O 9 M 
(Cia rdi et al.l I2OO0I ). However, the dim but abundant low 
mass sources (M < 1O 9 M ) may still str ongly influence 
the way in which the ionised reg ions grow ([Sokasian et al.l 
2003: IChoudhurv &; Ferrarall2007l ). These small galaxies are 
preferentially found in low-density environments along the 
cosmic web and may contain many metal- free Pop-Ill stars 
in the early Universe. It has been estimated that these low- 
mass halos account for about 80% of th e total ionising power 
at z ~ 7 ([Choudhurv &; Ferraral I2007P . Disregarding them 
may lead to an overestimate of the number of ph otons re- 
quired to reionise the Universe (| Sokasian et al.l 2003). 

Numerous theoretical studies have begun to investigate 
in detail the characteristic scales and the topology of the 
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cost and complexity of the radiative tran sfer problem, most 
simu la tions, with v e ry few exceptions (IGnedin &; Ost riker 
19971 : iKohler et all 120071 : IShin et al.1 120081 : IWise fe Abel 
2008), have treated reionisation through post-processing ap- 
plied to static or separately evolved gas density fields. This 
neglects the fact that the radiation field may exert important 
feedback effects on galaxy formation itself (e.g. llliev et al. 
120051 : lYoshida et aDl2007l : ICroft fe Altavll2008l ). It is there- 
fore an important task to develop more accurate theoretical 
models based on self- consistent simulations where the ion- 
isation field is evolved simultaneously with the growth of 
cosmic structures, a topic that we address in this work. 

In the literature, predictions based on simulations for 
the onset of reionis a tion range from re dshifts z ~ 30 — 
40 (llliev et al.1 120071: IWise fe Abel l2008h to z ~ 15 - 20 
(|Norman et al.1 1 19981 : lAbel et al.1 I2002T ). Reionisation then 
probably proceeds in an inhomo geneous and patchy fashion 
(iLidz et al.ll2007l : llliev et al J [20071 ), reflecting the inhomoge- 
neous density distribution of the large-scale structure. At 
first, many isolated ionised bubbles are formed. They then 
grow in size from ~ 1 Mpc during the earl y stages of reion- 



isation up to > 10 Mp c in the late phases (iFurlanetto et al 
20061 : llliev et al.1 [200(f). Around re dshift z ~ 13 (llliev et al 
2006), 8 < z < 10 (|Lee et aHl2008h , or z ~ 6 (|Gnedin fe Fan 
2006), the ionised regions overlap. These values are very 
uncertain and depend strongly on the modelling details of 
reionisation and the parameters of the underlying galaxy 
formation simulations. However, it is plausible that the fu- 
ture use of more self-consistent simulation techniques should 
be able to reduce the systematic modelling uncertainties. 

In this work, we therefore use the new radiative 
transfer algorithm we d eveloped in a previous study 
(|Petkova fe Springell2009h to carry out high-resolution sim- 
ulations of cosmic structure growth in the proper cosmo- 
logical context. The approximation to radiative transfer 
employed, the optically- thin variable Eddington tensor ap- 



proach ([Gnedin &; AbelfeoOlh , is fast enough to allow cou- 
pled radiative-hydrodynamic simulations of the galaxy for- 
mation process. At the same time, the employed moment- 
based approximation to the radiative transfer problem can 
be expected to be still reasonably accurate for the reionisa- 
tion problem. In particular, thanks to the photon-conserving 
character of our implementation of radiative transfer and of 
the chemical network, ionisation fronts are bound to prop- 
agate with the right speed. The Lagrangian smoothed par- 
ticle approach (SPH) we use automatically adapts to the 
large dynamic range in density developing in the galaxy for- 
mation problem. Combined with the fully adaptive gravi- 
tational force solver implemented in GADGET, this yields 
a numerical scheme that is particularly well suited for the 
cosmic structure formation problem. 

For a first assessment of our new approach we study 
simulations of the standard ACDM cosmology and treat 
star formation and supernova fee dback with the ISM 
sub-res olution model developed in ISpringel &; Hernquistl 
(2003a). For simplicity, we shall here only consider ordinary 
star- format ion regions as sources of ionising radiation. We 
are especially interested in whether the star formation pre- 
dicted by the simulations results in a plausible reionisation 
history of the Universe, and whether it at the same time 
yields a thermal and ionisation state of the IGM at interme- 
diate redshifts that is consistent with that probed by obser- 
vations of the Lyman-a forest. Finally, we are interested in 
possible differences induced in galaxy formation due to the 
spatially varying radiative feedback in the radiative transfer 
simulations, especially in comparison with the much sim- 
pler and so far widely adopted treatment where a spatially 
homogeneous UV background is externally imposed. 

This paper is structured as follows. We start in Section [2] 
with a brief summary our methods for simulating hydrogen 
reionisation. We then present our results in Section [3j focus- 
ing on the history of reionisation in Section l3TT1 the Lyman-a 
forest in Section 13.21 and the feedback from reionisation in 
Section 13.31 We end with a discussion and our conclusions 
in Section [4] 



2 SIMULATING HYDROGEN REIONISATION 

2.1 Radiative transfer modelling 

For our work we use an updat ed version of the cosmologi- 
cal si mulation code GADGET ([Springel et al.ll200ll ; ISpringel 
2005), co mbined with the radi a tive tr ansfer (RT) implemen- 
tation of IPetkova Springel ([2009). The RT equation is 
solved using a moment-based approach similar to the one 
proposed by IGnedin &; Abel (|200lh . The resulting partial 
differential equation essentially describes an anisotropic dif- 
fusion of the photon density field n 7 



<9n 7 

~dt 



1 dn^h 13 
k dxi 



(1) 



1 Differently from Petkova & Springel (2009) we solve equation 
fT|) for the photon density n 7 , rather than the photon over- 
density, with respect to the hydrogen density h 1 = n 7 /nn- 
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Figure 1. Ionised fraction in a slice through the simulated volume 
at evolution tim e t = 0.05 Myr of the cosmological density field 
test described in Iliev et al.l |20 



1 0° I" 1 1 1 I 



10 _1 - 



10" 



10~ 5 - 



10" 4 



10~ 5 - 



10' 



10" 7 




128 3 77=0.1 e = 6.4eV 
128 3 77=0.2 e = 6.4eV- 
128 3 t?=0.3 e = 6.4eV 
128 3 77=0.5 € = 6.4eV 
128 3 77=1.0 e = 6.4eV 



10 



12 



14 



redshift 



where c is the speed of light, k is the absorption coeffi- 
cient, h 13 is the Eddington tensor and s 1 is the source func- 
tion. The closure relation for this particular moment-based 
method is obtained by approximating the Eddington tensor 

h ij as 



h tJ = 



(2) 



(3) 



Tr(P) ' 

where P %J is the radiation pressure tensor 

ryij / \ / j3 / / /\ ( x — x )*( x — x )j 

P 3 (x)oc J d x p»(x ) _ — . 

This estimate of the Eddington tensors is carried out in the 
optically thin regime, giving the method its name. 

The abundances of hydrogen and helium species are up- 
dated by accounting for the processes of photo-ionisation, 
collisional ionisation and recombination: 

dnmi ~ ~ . ~ ~ ~ ~ , A , 

■ Jmnmn e nu + caonmn^nn - aHim e nminH, (4) 



dt 



dt 



dflHe 



= 7HemHem e riH 



7HeimHeim e ^H 
QiHelll^Helime^H, 



dt 



■ 7HeIinHeIin e riH - ttHeIimHeIim e nH, 



(5) 
(6) 

(7) 



where all abundances of hydrogen species (nm, ^hii) are ex- 
pressed in dimensionless form relative to the total number 
density nu of hydrogen, all helium abundances are fractions 
with respect to the helium number density, and the elec- 
tron abundance is expressed as h e = n e /nn. Furthermore, 7 
denotes the collisional ionisation coefficient, a is the recom- 
bination coefficient and <Jo = 6.3 x 10 -18 cm 2 is the photo- 
ionisation cross-section for hydrogen at 13.6 eV. These rate 
equations are solved using a semi-implicit scheme. 

The photo-heating of the gas by the ionising part of the 
spectrum is described by the heating rate 



nm 



-<j v (hv — hvo) 



nm cn^ em crm. 



(8) 
(9) 



Figure 2. Volume averaged neutral fraction as a function of 
redshift for the low resolution simulations. The different colours 
represent simulation results for different values of the efficiency 
parameter (or 'ISM escape fraction') 77 = 0.1, 0.2, 0.3, 0.5, 1.0. 
Reionisation happens earlier for higher efficiency since then more 
photons become available for ionising the gas. In all cases, the 
final phase of reionisation proceeds rapidly; over a small range of 
redshift, the neutral volume fraction drops from 10% to negligibly 
small values. 



where 



av— — <j v x 



4ttJ> 



(10) 



hv \ J hv 

is a frequency averaged photo-ionisation cross-section and 

4:7vl u a, 



dv^^^-(hv — hvo) x 
hv 



dv- 



hv 



(ii) 



is a frequency averaged photon excess energy (|Spitzerl ll998). 
For our calculations, we usually assume a black body spec- 



trum with T e ff = 10 K, which leads to am 



1.63 x 



10 -18 cm 2 and em = 6.4 eV. However, given the uncertain- 
ties in accurately calculating the net energy input during the 
photo-ionisation transition, which involves non-equilibrium 
processes, we will later on vary em in order to investigate 
the influence of different heating efficiencies. 

We have tested our radiative transfer implementation 
on the static cosmological density field that w as used in 
the ra diative transfer code comparison study bv llliev et al.l 
(2006). In Figure [1] we show the ionised fraction in a 
slice through the simulated volume at evolution time t = 
0.05 Myr. Reassuringly, our result is in good agreement with 
the ones obtained by other radiative codes in the comparison 
study. 



2.2 Treatment of star formation and radiative 
cooling 

We follow the radiative cooling of helium and hydrogen 
with a non-equilibrium treatment, with the rates as cited 
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Figure 3. Ionising background for the low resolution simu- 
lation as a function of redshift for different efficiency 77 = 
0.1, 0.2, 0.3, 0.5, 1.0. As expected, the background becomes 
higher for a larger efficiency since more of the photons emitted 
by the sources are made available to build up the ionising back- 
ground. The thin solid line shows the background computed by 
Haa rdt Madaul |l996) for quasars, which provides for an in- 
teresting comparison. Our results for the time evolution of the 
neutr al volume fr action agree quite well with previous studies 
(e.g.|Gnedin|[200a)- 

in IPetkova &; Springell ([20091 ). Winds and star formation in 
the dense, cold gas is modelled in a sub-resolution fash- 
ion (jSpringel fe Hernauistl l2003aJ ), where the interstellar 
medium is pictured as being composed of cold clouds that 
are embedded at pressure equilibrium in a hot tenuous phase 
that is heated by supernova explosions. Through the evap- 
oration of clouds, this establishes a tight self-regulation of 
the star formation rate. Previous work has shown that this 
model converges well with numerical resolution and yields 
star form ation rates that are consistent with the Kennicutt 
relation (|Kennicuttlll998h observed at low redshift. We shall 
here assume that the same star formation law also holds at 
high redshift. 

We use the star-forming regions in all simulated galax- 
ies as sources of ionising photons in our radiative transfer 
model. We adopt an ionising source lum inosity of iVsFR = 
10 53 photons M^yr (|Madau et al.lll999h . which relates the 
number of emitted photons to the star formation rate in 
units of M©yr _1 . This source luminosity is released individ- 
ually by every star- forming gas particle, hence the number 
of numerically represented sources is a non-negligible frac- 
tion of all simulation particles. Fortunately, the speed of our 
radiative transfer algorithm is almost insensitive to the total 
number of sources, because the Eddington tensor calculation 
can be carried out with a tree algorithm similar to the grav- 
ity calculation. This insensitivity of the computational cost 
to the number of sources is a significant advantage of the 
method used here, and is not shared by most alternatives 
for treating radiative transfer. 



To account for the uncertain absorption that occurs in 
reality in the spatially unresolved multi-phase structure of 
our simulations, we impose a phenomenological efficiency 
factor 77 on the source luminosities. In our simulation set we 
explore values in the range 77 = 0.1 — 1.0 to get a feeling 
for the sensitivity of our results to this uncertainty. We note 
however that 77 should not be confused with what is usually 
called galaxy escape fraction, which has a slightly different 
meaning. Our 77 is meant to be just an 'interstellar medium 
escape fraction' whereas photon losses in the gaseous halos 
of galaxies will be taken into account self-consistently in our 
simulations. 

2.3 Simulation set 

All our simulations assume a ACDM universe with cosmo- 
logical parameters Q = 0.27, Q\ = 0.73, Q b = 0.044, 
h = 0.7 and as = 0.9. In order to have sufficiently high 
mass resolution, we follow a comparatively small region in a 
periodic box of comoving size Lbox = 10/i -1 Mpc on a side. 
This region is still sufficiently large to give a representative 
account of the Lyman-a forest at redshift z = 3, which is 
the final time of our runs. 

Our primary simulation set has 2x 128 3 dark matter and 
gas particles, giving a mass resolution of 3.04 x 1O 7 /i _1 M 
and 5.29 x 10 6 /i _1 Mq in dark matter and gas, respec- 
tively. A selected model was also carried out at the much 
higher resolution of 2 x 256 3 , giving a gas mass resolution of 
6.62 x 1O 5 /i -1 M0. The gravitational softening was chosen 
as 1/35 of the mean particle spacing in each case, corre- 
sponding to e = 2.23/i -1 kpc and e = 1.12/i _1 kpc in the 
two resolution sets. In order to save computational time, 
most runs were restarted from z = 20 with different ra- 
diative transfer treatments, such that the higher redshift 
evolution did not have to be repeated. We have also sys- 
tematically varied the 'escape fraction' 77 and the heating 
efficiency e, considering the values 77 = 0.1, 0.2, 0.3, 0.5, 1.0 
and e — 6.4 eV, 16 eV, 20 eV, 30 eV for the low resolution 
set. For the high resolution run we chose the parameters 
77 = 0.2 and e = 30 eV. In this way we could systemati- 
cally determine the settings that give the most promising 
agreement with the observations. 



3 RESULTS 

3.1 Hydrogen reionisation history 

When star formation starts at around redshift z — 20 in 
our simulations, the process of reionisation begins through 
photo-ionisation of the gas around the star- format ion sites. 
However, as the dense gas has a high recombination rate, 
the progress in the reionisation is sensitively determined by 
a competition between the luminosity of the sources, the rate 
with which they turn on, and the density of the neutral gas 
they are embedded in. Using our simulation set, we first in- 
vestigate the global reionisation history and its dependence 
on the source efficiency parameter 77. 

Figure [2] shows the reionisation histories of our low res- 
olution simulated box for several choices of 77. As expected, 
the results for the reionisation redshift depend strongly on 
77. For larger efficiencies, the Universe gets reionised earlier, 
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Figure 4. Slices through the middle of the simulated volume for the low resolution simulation realisation with 77 = 0.1 and e = 30 eV. 
The maps are showing the density (left) with contours at ionised fraction nun = 0.001, 0.5 and 0.9 overlaid, ionised fraction (middle), 
and temperature (right) of the gas. The snapshots show a time evolution from top to bottom, with individual redshifts z = 7.2, z = 6.4, 
z = 5.7, z = 5 and z = 4, respectively. 
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Figure 5. Scatter plot of the neutral fraction as a function of over-density for two different redshift of the low resolution simulation 
realisation with 77 = 0.1 and e = 30 eV. The thick dashed line is the median neutral fraction. The neutral fraction shows a clear dependence 
with density, where high density regions are more ionised than low density regions. Very high density gas is highly ionised as well, because 
it is in the star forming stage. 




Figure 6. Lyman-a flux probability (left) and power spectrum (right) from the low r e soluti on si mulation w i th 77 = 0.2 and e = 6.4 eV at 
redshift z = 3. The result is compared with observational data from McDo nald et all ([20001 ) and lKim et a l. (2007). The flux probability 
agrees reasonably well with observations, whereas the power spectrum deviates at the high k end. We discuss this problem in the text. 



since more photons become available to ionise the hydrogen. 
This effect always overwhelms the reduction in star forma- 
tion and hence source luminosity that an increased efficiency 
parameter 77 also induces. In all cases the last phase of the 
Epoch of Reionisation (EoR) is rather short, i.e. the final 
10% of the volume transition very rapidly from being neu- 
tral to being essentially fully ionised. The redshift of reioni- 
sation when this occurs varies between z ~ 5 and z ~ 8 for 
77 = 0.1 and 77 = 1.0, respectively. In contrast, the EoR starts 
in general at much higher redshift. For example, the highest 
efficiency 77 = 1.0 model has already ionised nearly 30% of 
the volume by redshift z = 10. Interestingly, there is a sys- 
tematic variation of the time it takes to complete the last 



phase of the EoR. This is much more rapid for the high effi- 
ciency model than for the lower efficiency ones. Since typical 
bubble sizes at redshift z = 6 are up to 10/i _1 Mpc, we note 
that the size of our simulation box is really too small to draw 
any definitive conclusions about the global dura tion of reion- 
isation in the Universe. A number of authors (ICiard i et al. 
120031 : iFurlanetto fe Ohll2005l : iFurlanetto fe Mesingerl 12009) 
have pointed out that the local reionisation history depends 
on the environment, so that, for example, the evolution of 
the neutral fraction in a field or void region is different than 
in a proto-cluster region. As a result, only very large sim- 
ulation volumes of 100/i _1 Mpc or more on a side can be 
expected to yield a truly representative account of cosmic 
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Figure 7. Lyman-a flux probability (left) and power spectrum (right) for the low resolution simulation with efficiency 77 = 0.1, for four 
different va lues of the averag e d exce ss p hoton energy e = 6 .4 eV, 16 eV, 20 eV, 30 eV. The simulated data is compared to the observational 
result from lMcDonald et al.l (|2000h and lKim et"alT(|2007h . 
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Figure 8. Volume averaged neutral fraction as a function of red- 
shift. The set of low resolution simulations with efficiency 77 = 0.1 
and average photon excess energy e = 6.4 eV, 16 eV, 20 eV, 30 eV 
is compared with the high resolution run with 77 = 0.2 and 
e = 30 eV. When a higher energy per ionising event is injected, 
the universe gets ionised slightly earlier since higher temperatures 
help to maintain higher ionised fractions. 

reionization. Since at the same time an equal or better mass 
resolution as we use here is required, such calculations are 
very expensive and beyond the scope of this work, but the 
should become feasible in the future on the next generation 
of high-performance computers. 

In Figure [3] we compare the ionising background for the 
same simulations, which gives further interesting clues about 
the reionisation histories of the different models. The back- 
ground is computed as the volume- averaged intensity in the 
simulation box. The ionisin g background is c ompa red with 
an analytical estimate by iHaardt &; Madaul (| 19961 ) for the 
meta-galactic ionising flux from quasars rather than stellar 
sources, which is an interesting comparison point for the 
expected value of the background. Clearly, lowering the ef- 
ficiency 77 decreases the ionising background as well, con- 
sistent with the findings above. Interestingly, the rapid rise 
of the mean background intensity ends when reionisation is 
complete. From this point on, the background shows only a 
weak residual evolution. 

The time evolution of the temperature, ionised frac- 
tion and density fields in a representative simulation model 
(77 = 0.1 and e = 30 eV) is illustrated in Figure H The dif- 
ferent panels correspond to slices through the middle of the 
simulation box, between redshifts z — 7.2 and z — 4, from 
the top row to the bottom row. It is seen that the ionised 
regions start to grow first around high density peaks, where 
the star forming regions are concentrated. Then the radi- 
ation diffuses into the inter-cluster medium. The filaments 
remain less ionised than the voids for a while, since their 
density is much higher. Initially the photons heat up the gas 
in the ionised region to temperatures slightly above 10 4 K. 



Figure 9. Ionising background as a function of redshift. The set 
of low resolution simulations with efficiency 77 = 0.1 and averaged 
photon excess energy e = 6.4 eV, 16 eV, 20 eV, 30 eV is compared 
with the high resolution run with 77 = 0.2 and e = 30 eV. For a 
higher injected energy per ionisation event, the background also 
increases as a result of the higher temperature, which leaves more 
photons unabsorbed so that they can contribute to a higher level 
of the ionising background. 

As the ionising background declines due to the expansion 
of the Universe and the drop of star formation, the heat- 
ing becomes less effective and the temperature of the highly 
ionised gas in the voids drops somewhat as a result of the 
expansion cooling. 

In Figure [5j we present a contour plot of the neutral 
fraction versus the over-density of the gas in the represen- 
tative model for two different redshifts, at z — 6.4 before 
reionisation is completed, and at z — 3 after reionisation 
is completed. There is a clear dependence of the neutral 
fraction on over-density in both cases. The high density re- 
gions around star-forming matter are ionised very quickly. 
The average density regions, e.g. filaments, tend to be less 
ionised and get ionised after the low density regions. How- 
ever, note that in the star-forming tale of the diagram all 
the gas is ionised. This is here due to the star formation 
scheme adopted in GADGET, where star- forming gas parti- 
cles are assigned a mean mass-weighted temperature which 
is so high that all this gas is formally collisionally ionised. 
We note t hat these res ults are very similar to the ones re- 
ported bv iGnedinl (2000), where a similar relation between 
neutral fraction and density is found. However, we are able 
to probe somewhat higher densities thanks to better spatial 
resolution of our simulations. 

3.2 Lyman-a forest 

We next turn to an analysis of the thermal state of the inter- 
galactic medium left behind at z — 3 by our self- consistent 
reionisation simulations. To this end we use the gas den- 
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sity, gas temperature, gas velocity and ionisation state of 
the gas in the simulation box and compute Lyman-a absorp- 
tion spectra for random lines of sight. We then compare the 
statistics of these artificial absorption spectra with observa- 
tional data on th e Lym an- a forest at th i s epoc h, as given by 
iMcDonald et al.1 (|2000h and iKim etHI (|200<f ). Here we as- 
sume that the main source of ionising sources is stellar and 
thus discard any c ontribution from a quasar- type spectra, 
in agreement with iMadau et al.l (|l999h . We also note that 
some discrepancies are possible due to the fact that helium 
is not photo-ionised in our simulations, but only collisionally 
ionised. 

Our simulations have the necessary gas mass resolution 
at redshift z = 3 required to reprodu ce realistic Lyman-o; ab - 
sorption in the low density regions ([Bolton Bec ker 2009). 
However, we can not match the required simulation volume 
of size ~ 40/i _1 Mpc to properly sample the largest voids. 
This can have a significant effect on our predictions of the 
Lyman-a flux probability distribution and power spectra. 

Figure [6] shows the flux probability distribution function 
(PDF) and flux power spectrum for our simulation model, 
where an efficiency parameter of r/ = 0.2 and averaged pho- 
ton excess energy e = 6.4 eV were adopted. For this choice, 
we achieve the best fit to the flux PDF. However, for all 
the other models the power spectrum is over-predicted at 
high wave numbers. We suggest that this overestimation of 
the power spectrum is due to the insufficient heating of the 
gas in low density regions, causing an excess of small-scale 
structure in the Lyman-a forest. 

In order to examine this effect further, we vary the pho- 
ton excess energy e used in the photo-heating and examine 
the influence this has on the flux probability and power spec- 
trum. There are two possible reasons why our simulations 
underestimate the photo heating. First, we expect that some 
non-equilibrium effects in the photo-heating are treated in- 
accurately due to our implici t treatment of the rad iation 
transport and chemistry (e.g. iBolton fc Beck er 2009). sec- 
ond, photo-heating is different in optically thin and opti- 
cally thick regions. For example, in an optically thick region 
the average photon excess energy obtained from Eqn. 
is e = 29.9 eV. It is however likely that our approximative 
radiative transfer scheme leads to inaccurcies in the effec- 
tive heating rates of regions of different optical depths, due 
to the varying accuracy of the scheme in different regimes. 
Part of these inaccuracies can be absorbed into a suitably 
modified value of the effective heating rate e. To explore the 
full range of plausible values, we therefore vary the values for 
e as follows: e = 6.4 eV, 16 eV, 20 eV and 30 eV. We aim to 
bracket what can be expected when non-equilibrium effects 
are fully taken into account in future treatments, and want 
to identify the case that provides the best representation of 
the Universe at redshift z = 3. 

Figure [6] shows the flux PDF and power spectra for 
these different heating values. Clearly, the high wave num- 
ber region of the flux power spectrum is strongly influenced 
by the amount of injected heat energy into the gas, and the 
increase of the temperature also affects the flux probabil- 
ity distribution. For the low efficiency of 77 = 0.1, there 
is a substantial mismatch already in the flux PDF, sim- 
ply because there is too little ionisation overall so that the 
mean transmission ends up being too low. However, as the 
adopted photo-heating energy increases, the gas is getting 
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Figure 10. Star formation rate density as a function of redshift 
for the low resolution simulation set at different efficiencies of 
77 = 0.1, 0.2, 0.3, 0.5, 1.0. The results are compared to the SFR 
history of a low resolution simulation with in stantaneous reioni- 
sation at z = 6 and photo- heating by a lHaardt fc M adau dl996h 
ionising background (thin black line), and a simulation with nei- 
ther reionisation nor photo-heating (thick black line). The photo- 
heati ng from stellar sources decreases star formation, as suggested 
bv lPawlik fcSc have (2009). As the escape efficiency gets higher, 
this effect becomes progressively stronger. 



hotter and is able to stay ionised longer due to the higher 
temperatures, yielding a better fit to the flux PDF. At the 
same time, small-scale structure in the flux power spectrum 
is erased due to thermal broadening, bringing the simula- 
tions into agreement with the observation. This shows the 
power of detailed Lyman-a data to constrain simulations of 
the reionisation process. In our current models we need to 
adopt a quite extreme heating efficiency of 30 eV combined 
with a low 'escape fraction' of 77 = 0.1 to achieve a good 
match to the data. 

In Figures [8] and [9] we compare the impact of the differ- 
ent photo-heating efficiencies on the evolution of the neutral 
volume fraction and the ionising background. We also show 
for comparison the results from our high resolution simu- 
lation, which is discussed below in the text. As expected, 
an increase in the heating energy leads to a slightly earlier 
reionisation and to a slightly elevated ionising background 
flux. Both of these effects can be readily understood from 
the higher gas temperature produced in the ionised gas when 
the higher heating efficiency is adopted. However, the effect 
is quite weak, and very much smaller than the changes re- 
sulting from a different choice of 77. 

We have also measured the Thomson electron scattering 
optical depth in our high resolution simulation and found it 
to be T es = 0.049, whic h is smaller than the WMAP7 value 
t wmap = .088±0.015 (|Komatsu et alJl2010h . This discrep- 
ancy is, however, not critical since the simulated volume is 
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Figure 12. Lyman-a flux probability (left) and power spectrum (right) for the high resolution simulat ion w ith efficiency 77 = 0.2 and 
averaged excess photon energy e = 30 eV, compared to observational results from lMcDonald et al.l (|2000T ) and lKim et al.l (|2007l ) . 
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Figure 11. Star formation rate density as a function of redshift. 
The set of low resolution simulations with efficiency 77 = 0.1 and 
averaged photon excess energy e = 6.4 eV, 16 eV, 20 eV, 30 eV are 
compared to the high resolution run with 77 = 0.2 and e = 30 eV. 
The results are compared to the SFR history of a low resolution 
simulation with instantaneous reionisation at z = 6 and photo- 
heating by a lHaardt fc Madaul Jl996|) ionising background (thin 
black line), and a simulation with neither reionisation nor photo- 
heating (thick black line). The star formation decreases with in- 
creasing heating energy, as expected. For the low resolution run 
with 30 eV, the result of the self-consistent radiative transfer cal- 
culation matches the simulation with instantaneous reionisation. 
The high resolution simulation SFR is higher at higher redshift 
due to better resolution and agrees well with the other results at 
redshifts less than z = 6. 



too small to obtain a realistic value and we have also not 
included photo-ionisation of helium. 



3.3 Feedback from reionisation 

Reionisation can in principle exert a strong feedback effect 
on the gas through the temperature increase induced by 
photo-heating. As the gas temperature increases, the gas 
densities will be lowered through pressure effects. The gas 
will then cool and collapse more slowly, such that the star 
formation rate is ultimately reduced. Especially small dark 
matter halos should be sensitive to this effect. In the extreme 
case of halos that have virial temperatures comparable to or 
only slightly larger than the temperature reached by the gas 
through reionisation, the UV radiation may even completely 
suppress atomic cooling and efficient star formation. This 
effect is often invoked to explain why so many of the dark 
matter satellites expected in ACDM in the halos of ordinary 
L* galaxies are apparently largely devoid of stars. 

In order to highlight the radiative transfer effects on the 
star formation in galaxies, we compare our simulations with 
two fiducial models where no radiative transfer is used. The 
first is a simulation simply without any photo-heating of 
the gas, while the second one is a simulation where reionisa- 
tion is induced by an externally imposed, spatially homoge- 
neou s UV background based on a modified Ha ardt &; Madaul 
(1996) model that causes reionisation and an associated 
photo - heating of the gas at z = 6 (for details see lDave et al.l 
1999). The latter model corresponds to the standard ap- 
proach applied in many previous h ydrodynamical simula- 
tion models of galaxy forma tion (e.g. iTornatore et alJ l2003: 
IWadeouhl fe SpringelllioToh . 

In Figure [TOl we compare our results for the cosmic star 
formation rate density evolution as a function of the adopted 
efficiency parameter 77 (for the low resolution simulation set). 
We also include the two fiducial comparison models as lim- 
iting cases. As we increase the escape efficiency of the ion- 
ising radiation, the star formation drops, as expected, since 
this makes more photons available to photo-heat the gas. 
We note that our results for the SFR are always lower than 
the fiducial simulation where no photo-heating is i n clude d 
at all, consistent with findings bv IPawlik &; Schavel (j20Q9h . 
Towards lower redshift, the reduction of the SFR due to the 
radiation field becomes progressively larger. The run with 
77 = 0.2 quite closely corresponds to the simulation with the 
imposed reionisation epoch, but starts to slightly differ at 
redshifts z < 4. 



© 2010 RAS, MNRAS 000, HHH] 



Simulations of galaxy formation with radiative transfer 11 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 



10 4 



10" 



10° ■ 



10" 




log(p/</)>) = -1, without photoheoting" 



log(p/</0>) = - 1, with photoheotfng 
Iog(p/<y0>)=0, without photoheoting 
log(p/<p>)=0, with photoheoting 
log(p/</0>)=1, without photoheoting 
log(/)/<y0>)=1, with photoheoting 

■ ■ i i i i i ■ 

3 4 5 6 7 

redshift 




log(1 + 6) 



Figure 13. Evolution of the mass averaged temperature at three 
different over-densities log 10 (p/(p)) = —1,0,1 for the low reso- 
lution simulation with 77 = 0.1 and e = 30 eV, with and without 
photo-heating. The strongest effect is observed in the low density 
gas, which is heated by the photons much more than the higher 
density gas. At all densities, however, photo-heating increases the 
temperature, as expected. 



We also carry out a corresponding comparison for a low 
resolution simulation set with constant efficiency 77 = 0.1 
but different values for the photon excess energy. In Fig- 
ure [TT] we show the results for the SFR, again including the 
two fiducial models as limiting cases for comparison. The 
results confirm the expectation that an increase of the pho- 
ton excess energy decreases the star formation rate density. 
Interestingly, the model that best reproduced the Lyman-a 
power spectrum observations, the one with e = 30 eV, quite 
closely follows the star formation rate density obtained for 
the fiducial model where reionisation is imposed at z = 6. 

For our high resolution run we chose to repeat the sim- 
ulation with averaged photon excess energy e — 30 eV and 
adopt a higher escape fraction 77 = 0.2 rather than 77 = 0.1. 
In this way we make sure we account for the trapping of 
photons in high density peaks, which were not present in 
the low resolution runs. In Figure [TT] we show how the star 
formation rate history compares to the ones from the low 
resolution runs. They are in good agreement, except for the 
higher redshift, where the high resolution captures mor e 
star formation, as expected (|Springel &; Hernquistll2003bh . 
The Lyman-a forest flux probability and power spectrum at 
z — 3 for this simulation are shown in Figure 1121 While the 
simulated data is in reasonable qualitative agre e ment with 
the observationa l results from iMcDonald et al.l (|2000) and 
iKim et all (|2007h , it does not provide in this case a detailed 
fit within the error bars, again highlighting that simultane- 
ously accounting for the cosmic star formation history, cos- 
mic reionisation and the state of the IGM at intermediate 
redshifts provides a powerful constraint on self- consistent 
simulations of galaxy formation and reionisation. 



Figure 14. Median temperature of the gas as a function of over- 
densities log 10 (l + 5) at redshifts z = 5.09 and z = 3 for the high 
resolution simulation with 77 = 0.2 and e = 30 eV. At redshift z = 
5.09, shortly after reionisation is completed, the temperature at 
low densities is clearly higher than that at higher densities (apart 
from the gas in the star-forming phase). This can be interpreted 
as an inverted equation of state. At lower redshift the relation 
reverts again to normal form as the gas in the low density regions 
cools down adiabatically due to the expansion of the Universe. 
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Figure 15. Mean stellar and gas masses as a function of the DM 
halo mass at z = 3 in the high resolution simulation. The black 
vertical corresponds to a mass of 100 DM particles, which can be 
taken as an (optimistic) resolution limit of the simulation. Photo- 
heating slows down the collapse of gas in halos, which in turn also 
decreases their stellar and gas masses. The effect becomes stronger 
for low mass DM halos. 
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Figure 16. Baryon fraction as a function of the DM halo mass 
at z = 3 in the higher resolution simulation. The black vertical 
corresponds to a mass of 100 DM particles. 
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Figure 17. Baryon fraction as a function of DM halo mass at 
z = 3 in the lower resolution simulations with efficiency 77 = 0.1, 
and for different averaged photon excess energy. The dashed line 
shows the baryon fraction when no photo-heating has taken place. 
The black vertical line corresponds to a mass of 100 DM particles. 
Photo-heating does not affect the baryon fraction as strongly as in 
the high resolution simulation. We also observe that differences in 
the excess photon energy do not have a large effect on the baryon 
fraction. 



In Figure [13] we explore the temperature evolution of 
the gas at different characteristic densities, corresponding to 
under- dense gas by a factor of 10, gas at the mean density, 
and over-dense gas by a factor of 10 relative to the mean. 
We compare our default simulation with radiative transfer 
and photo-heating to the fiducial simulation where no such 
heating is included at all. Clearly, the effect of photo-heating 
is most prominent in the lowest density gas. This gas is only 
weakly heated by structure formation shock waves when 
photo- heating is not included. In contrast, when reionisa- 
tion is accounted for, the temperature of this gas reaches 
a high value of ~ 10 4 K at the end of the epoch of cosmic 
reionisation, and even before that, the mean temperature of 
this gas is raised considerably as a result of the patchy and 
temporally extended reionisation transition in our radiative 
transfer simulations. Interestingly, after reionisation is com- 
plete, the mean temperature of the under-dense gas starts 
to slowly decline again, while already for the mean density 
gas structure formation shocks can provide for a slow further 
increase of the temperature. 

We also analysed the median temperature of the gas 
as a function of over-density. As shown in Figure 1141 after 
reionisation has been completed, the low density gas ends 
up with a higher median temperature than the higher den- 
sity gas (except for the gas in the star- forming phase) . This 
points towards an 'inverted equation o f sta te', as observed by 
iBolton et all (|2008T ). lTrac et al.1 (|20Q8h and lFurlanetto fe Ohl 
( 200£j). At later times, the equation of state reverts again to 
a normal positive slope, when the low density gas cools down 
due to the adiabatic expansion of the Universe. 

Finally, in Figure 1151 Figure [16] and Figure [17] we ex- 
plore the impact of the ionising radiation field on the gas 
and stellar mass content of individual dark matter halos. 
To this end we run a group finder on our simulations and 
simply determine the average gas mass, stellar mass and 
baryon fraction of halos as a function of their dark matter 
mass. We compare the z = 3 results of our higher resolu- 
tion radiative transfer simulation with the simulation where 
photo- heating is completely ignored. Interestingly, we find 
a reduction of the gas and stellar mass for all halo masses 
when radiative transfer is included. The effect is quite weak 
for large halos but becomes progressively larger for small ha- 
los. At dark matter halo masses of Mdm = 10 9 M the sup- 
pression in baryonic content is approximately 60%, while at 
Mdm = 10 12 M0 it drops to only a few percent. This shows 
clearly the important impact of the ionising radiation field 
on small dwarf galaxies, in particular. While an externally 
imposed UV background can perhaps acc ount for the mean 
effect of th is radiative feedback process (Hoe ft et al.i r2006: 
lOkamoto et al.ll2008h , only a spatially resolved treatment of 
radiative transfer can account for effects of proximity that 
may well play an impor tant role in shapi n g, e.g., the satellite 
luminosity func tion (|Munoz et alJl2009l : Busha et alJl20ld : 
llliev et alJl2O10h . 



4 DISCUSSION AND CONCLUSIONS 



We have p resented the first applicat ion of our new imple- 
mentation ([Petkova &; Sprin gel 2009) of radiation hydrody- 
namics in the cosmological simulation code GADGET. We 
focused on the problem of cosmic reionisation, aiming in 
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particular at a first test on whether the default star forma- 
tion model in the code combined with our radiative transfer 
modelling can yield a plausible reionisation history of the 
Universe and a reasonable thermal state of the intermedi- 
ate redshift intergalactic medium. For simplicity, we have 
here only studied star-forming galaxies as ionising sources, 
and restricted the analysis to hydrogen reionisation alone. 
Based on the encouraging results collected here, it is clearly 
worthwhile to extend the model further in future work. 

Since the level of internal absorption in the interstellar 
medium is uncertain and cannot be resolved by our sim- 
ulations, we have examined models with different effective 
source efficiencies r\. Likewise, as we have not included a 
detailed spectral treatment and the time evolution of non- 
equilibrium in the chemistry may be inaccurate, we have 
parametrised the heat input per ionisation event in terms of 
a parameter e. 

We find that our simulated universes can get reionised 
by star formation in ordinary galaxies alone, with the epoch 
of reionisation ending between redshifts z = 8 to z — 5, de- 
pending on the assumed escape efficiency. The final phase 
transition is always quite rapid in this our setup, but sizable 
fractions of the volume begin to be reionised much earlier. 
The heating efficiency has only a weak influence on the reion- 
isation history, but a stronger one on the cosmic star forma- 
tion rate density. In fact, we have shown that photo- heating 
plays an important role in the evolution of the baryonic gas. 
As a result of the associated heating, it changes baryonic 
structure formation by slowing down the collapse of gas in 
dark matter halos, thereby delaying star formation. This ef- 
fect is strongest for the lowest mass halos, where the DM 
potential well is not deep enough to easily overcome the 
thermal pressure from the effects of photo-ionisation. 

Our simple models of a self-consistent treatment of 
galaxy formation and radiative transfer are not only able 
to produce a plausible history of reionisation, but they also 
manage to approximately match the basic statistics of the 
Lyman-a forest, at least for an appropriate choice of the pa- 
rameters 77 and e. This suggests that the low-redshift IGM 
data can be a powerful additional constraint on future reion- 
isation modelling in structure formation simulations. 

Despite these encouraging results it is also clear that our 
simulation results are likely still affected by numerical reso- 
lution effects, because the resolution in the lowest mass halos 
is still too coarse to yield fully converged results. Ideally, we 
would like to resolve the full range of star-forming halos with 
enough particles to achieve fully converged results. While 
this is unlikely to qualitatively change any of the results 
presented here, future precision work will require such cal- 
culations. Another important caveat that will require further 
study are uncertainties due to the radiative transfer approx- 
imation itself. This is probably best addressed by comparing 
the results with a completely different approach to radiative 
transfer. We are presently finalising such a scheme in the 
moving- mesh code AREPO. It avoids the use of the diffusion 
approximation and can cast sharp shadows, hence a direct 
one-to-one comparison with the optically-thin variable Ed- 
dington tensor approach applied here will be particularly 
interesting. 
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